DISCONTINUOUS GALERKIN METHODS 
FOR MASS TRANSFER 
THROUGH SEMI-PERMEABLE MEMBRANES 
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Abstract. A discontinuous Galerkin (dG) method for the numerical solution of initial/boundary value multi- 
compartment partial differential equation (PDE) models, interconnected with interface conditions, is presented and 
analysed. The study of interface problems is motivated by models of mass transfer of solutes through semi-permeable 
membranes. More specifically, a model problem consisting of a system of semilinear parabolic advection-diffusion- 
reaction partial differential equations in each compartment, equipped with respective initial and boundary conditions, 
is considered. Nonlinear interface conditions modelling selective permeability, congestion and partial reflection are 
applied to the compartment interfaces. An interior penalty dG method is presented for this problem and it is analysed 
in the space-discrete setting. The a priori analysis shows that the method yields optimal a priori bounds, provided 
the exact solution is sufficiently smooth. Numerical experiments indicate agreement with the theoretical bounds and 
highlight the stability of the numerical method in the advection-dominated regime. 

1. Introduction. Models of mass transfer of substances (solutes) through semi-per- 
meable membranes appear in various contexts, such as biomedical and chemical engineering 
applications |2o*l . Examples include the modelling of solute dynamics across arterial walls 
(see, e.g., [52 1 and the references therein) and cellular signal transduction (see, e.g., (14 1 and 
the references therein). 

This work is concerned with the development and analysis of numerical methods for a 
class of continuum models for mass transfer based on initial/boundary value multi-compart- 
ment partial differential equation (PDE) problems, closed by nonlinear interface conditions. 
The interface conditions considered are the Kedem-Katchalsky (KK) equations, which repre- 
sent an established model for the mass transfer mechanisms ll37l [361 . More specifically, we 
consider a generic model problem consisting of a system of semilinear advection-diffusion- 
reaction parabolic PDE problems in multi-compartment configurations, coupled with nonlin- 
ear interface KK-type conditions. The focus is to address some challenges in the numerical 
solution of these models, such as the treatment of non-linearities due to both the interface 
modelling and the non-linear reactions, the discontinuity of the state variables across the in- 
terface, as well as the development of stable numerical methods in the advection-dominated 
regime. 

Numerical methods for mass transfer problems based on conforming finite elements have 
been developed for the solution of solute dynamics across arterial walls; see l52l |451 l44l 
and the references therein for more details. Some existence results for the purely diffusing 
interface problem without forcing, coupled with KK-type interface conditions, along with 
some numerical experiments are given in lfl3ll . Further, numerical approaches to the treatment 
of interface conditions for PDE problems, resulting to globally continuous solutions can be 
found, e.g., in @ EJ Q3] M SD • 

Discontinuous Galerkin (dG) methods (see, e.g., PT1 l46l [181 and the references therein) 
and mortar methods (see, e.g., (9) for a survey) have also been proposed for the treatment of 
coupled systems via interface conditions in various contexts [49 27, 28 l22l l23l [33l [30l l29l . 
Also, the advantages of dG methods for interfacing different numerical methods (numeri- 
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cal interfaces) have been identified ll42l[T7l . as well as their use on transmission-type/high- 
contrast problems, yielding continuous solutions across the transmission interface, has been 
investigated l2ll[TTll24l[T2l. 

Here, we consider a dG method of interior penalty type for the solution of the semilinear 
parabolic advection-diffusion-reaction PDE system coupled with nonlinear interface condi- 
tions of KK-type across the subdomains. The use of dG is motivated partly by the observa- 
tion that the interface conditions, yielding discontinuous solutions across the interface, can 
be imposed by modifying the interior penalty dG numerical fluxes. Another important factor 
for employing a dG method is the desired stability property of the numerical method in the 
advection-dominated regime. A priori bounds for the proposed spatially discrete dG method 
in both the i°°(I/ 2 )- and L 2 (7? 1 )-type norms are presented for a range of reaction fields, 
under the simplifying assumption that the finite element mesh is aligned with the subdomain 
interfaces. 

A priori error bounds for interior penalty dG methods for parabolic problems have been 
considered in various settings (see, e.g., Il46ll for an exposition and more recently [19]). DG 
methods for semilinear parabolic spatially self-adjoint problems with locally Lipschitz con- 
tinuous nonlinearity have been analysed in ll40l . In the present analysis, advection terms 
are included and systems of equations are considered. In the presence of advection, the 
analysis of the symmetric interior penalty dG method in [40 1 would require the assumption 
of quasi-uniformity of the mesh. To avoid this assumption, a different continuation argu- 
ment is employed in the derivation of the a priori bounds presented here, at the expense of a 
stricter growth condition on the nonlinearity of the forcing term. This continuation argument 
is inspired by the derivation of a posteriori bounds for semilinear parabolic phase-field mod- 
els (38 1 8 1. The non-linear interface terms are tackled using a non-standard elliptic projection 
which is inspired by a classical construction of Douglas and Dupont [20 1 for the treatment of 
non-linear boundary conditions. Further a priori error estimates under the weaker assumption 
of locally Lipschitz strength of the reaction term, at the expense of stricter mesh assump- 
tions, are also shown. The fixed point argument used has been applied to other types of finite 
element methods for time-dependent semilinear problems, cf. Il2l l25l . 

The paper is organized as follows. In Section [2] the PDE model along with a short 
derivation of the non-linear interface conditions is presented. Section [3] is devoted to the 
description of the dG method proposed for the advection-diffusion part of the spatial operator 
incorporating the non-linear interface conditions, while Section|4]contains error estimates for 
the new elliptic projection, which will, in turn, be utilized as the elliptic projection in the 
subsequent analysis. Two a priori error bounds are presented in Section [5] while Section [6] 
contains numerical experiments highlighting the stability and the optimal rate of convergence 
of the proposed method in practice. Finally, we draw some conclusions in Section [7] 

2. Interface modelling and governing PDEs. We shall consider systems of parabolic 
semi-linear PDEs describing the flux of solutes around and through a semi-permeable mem- 
brane. The membrane is modelled as an internal boundary equipped with nonlinear interface 
conditions which are described in the following section. 

2.1. Interface modelling. We outline the Kedem-Katchalsky (KK) equations modelling 
solutes flow across semi-permeable membranes. The KK equations have been introduced 
in [37]; we refer to [50| for earlier works and to [26 1 for a thorough exposition. 

It is assumed that the membrane separates two compartments filled with a free fluid 
which is called the solvent and that the membrane characteristics are uniform in space. The 
KK equations specify the dependence of the solutes and solvent fluxes across the membrane 
in terms of two driving forces, namely the hydrostatic and osmotic pressure jumps. In the 



DG METHODS FOR MASS TRANSFER THROUGH MEMBRANES 3 

case of a single solute, the solvent flux J v and the solute flux J s are given by 

J v = L P (Sp - a d 6ir), (2.1) 

J s = u) 5tt + J v (l — <Jf)u, (2.2) 



in terms of the hydrostatic and osmotic pressure jumps 5p and 5tt — RTSu (with R being the 
ideal gas constant and T denoting temperature). Here, u represents the average concentration 
of the solute across the membrane, and Lp, a d and at, and cj are the phenomenological 
membrane transport coefficients of filtration, reflection, and permeation, respectively. 

Equation \2.\\ , which is known as Starling's law of filtration, shows that the solvent flow 
is affected by the osmotic flow of the solute. This observation introduces a nonlinearity in the 
transport term of the solute flux equation \2.2\ . Indeed, substituting J v into ( |2.2[ ) we get the 
final model for the solute flux 

J s = u) RT 5u + Lp(5p — adRTSu)(l — af)u = p(u) Su + rbu. 

In this last expression, we have collected the diffusive part of the flux by introducing the 
non-linear permeability function p(u), and written the advective transport in terms of the 
friction coefficient r € [0, 1] and the normal component of the transport field yielded by the 
hydrostatic pressure, which is denoted by b. 

It is necessary to fix a model for the average concentrations inside the membrane. In 
the case of relatively thick membranes ll52l l44l . it is appropriate to consider a mechanical 
approach describing solutes flow within the membrane through an advection-diffusion model. 
This approach yields a weighted arithmetic average of the concentration u 1 and u 2 at the 
membrane's faces: 

u = 5 w u := w l u l + w 2 u 2 , (2.3) 

for some given weights w = (w , w 2 ) with w 1 + w 2 = 1. These can be expressed in terms 
of the internal Peclet number of the membrane advection-diffusion model and are so that the 
upwind value dominates [52], cf. the conditions given below in ( |2.12[ ). 

Usually, in the case of n solutes, it is assumed that there is no direct coupling between 
the solutes' fluxes ll26l . Under this simplifying assumption the KK equations become 

n 

J v = L P (5p - ^ a d,j Sttj) (2.4) 
j=i 

n 

J s ,i = Pijuj) Suj - / ]pi,j{uj) Suj +rjbuj Viel, ...,n, (2.5) 
j'=i 

with Ui representing the concentration of the i-th solute, pi the corresponding permeability, 
and pij the cross-coefficients expressing the crowding effect inside the membrane. In what 
follows, the fluxes given by ( |2.5| l are used to close the PDE problem with appropriate interface 
conditions. 

2.2. Notation. We denote by L p (ui), 1 < p < +oo, the standard Lebesgue spaces, 
uj C R rf , d E {2,3}, with corresponding norms || • \\ P;U '< if w = fi, we shall write instead 
|| • ||p. Also the norm of L 2 (u>) will be denoted by ||-|| w and if cj = £1 by ||-|| for brevity; 
by (•, •) we write the standard L 2 -inner product on f2; when the arguments are vectors of 
L 2 -functions, the L 2 -inner product is modified in the standard fashion. We denote by H s (cu) 
the standard Hilbertian Sobolev space of index s € M. of real-valued functions defined on 
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ft 1 



FIG. 2.1. The domain of solution Q is subdivided into two subdomains f2 , D. 2 . The interface boundary is 
defined asT 3 := (dU 1 ndCl 2 )\dn. 



uj C M. d ; in particular Hq(lj) signifies the space of functions in -ff 1 (w) whose traces onto the 
boundary dtu vanish. For 1 < p < +oo, we also define the spaces L p (0, T; X), with X being 
a Banach space with norm || ■ \\x, consisting of all measurable functions v : [0, T] —> X, for 
which 

\\v\\lp(o,t-,x) : = (J q IK*)Hx d *) ^ < for 1 < P < 

\\v\\l^(o.T;X) ■= css sup ||u(£)||x < +oo, for p = +oo. 

0<t<T 

Finally, we denote by C(0,T;X) the space of continuous functions v : [0, T] — > X with 
norm |M|c(o,t ; x) := max < t < T < +o°- 

Let f2 be a bounded open domain with Lipschitz boundary in R d , and let dfl be the 
boundary of ft. The domain ft is subdivided into two subdomains il 1 and fl 2 , such that 



fl = n 1 U ft 2 U r J; where Tj := {dVL 1 n dn 2 )\dn, see Figure \2A\ For i = 1,2, we 
assume that Q l has Lipschitz boundary and that dft 1 D dfl has positive (d — 1) -dimensional 
(Hausdorff) measure. We define H s := [H^fl 1 U n 2 )] n , sei 

We shall employ the following notational convention: vectors are indicated with lower 
case bold symbols, n x n diagonal matrices with upper case (non-bold) symbols, and n x d 
tensors with upper case bold symbols. 

The gradient Vv of a vector function v : SI 1 U fi 2 — > W 1 in H 1 is a mapping ft 1 U il 2 — > 
M. nxd gained from componentwise application of the gradient operation: 

Vv := (Vv 1 ,...,Vv n ) T . 

Similarly the divergence V • Q of the tensor-valued function Q : O 1 U S7 2 — > W ixd is 
V • Q := (V • Qi, . . . , V • Q n ) T where the Qi are rows of Q. 

2.3. Model Problem. For a time interval [0, T], T > 0, and for 

u e L 2 (0, T\ H 1 ), u t € L 2 (0, T; H _1 ), 

with u := (iti , . . . , u„) T , we consider the system of semilinear parabolic equations 

u ( -V-(AVu-[/B) + F(u) = in (0,T] x (fl 1 U ft 2 ), (2.6) 

with U denoting the diagonal matrix U = diag(ui, . . . , u n ). Here, B is an n x d tensor field 
with rows Bi £ C 1 (0,T;W^ 1 ' oo (ft\r J ) d nVK oo (div,ft)),i = 1, . . . , n, and A € [C*([0,T] x 
ft 1 U Q 2 )] nxn diagonal, with A = diag(a x , a 2) . . . , o n ), where a, : [0, T] x ft 1 U £! 2 — s- R, 
i = 1, . . . ,n. We assume that there exists a constant a m i n > of uniform parabolicity such 
that dj(t, x) > a m i n for alii = 1, . . . , n and (t, x) € [0, T] x ft. For simplicity, we also 
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require that the matrix 1/2 diag(V • B) is positive semi-definite. Finally, F : R" — > K." is a 
vector field satisfying the growth condition 

|F(w) - F(v)| < (7(1 + |w| + |v|) 7 |w - v|, (2.7) 

for w, v G R n and 7 > constant, where | • | denotes the Euclidean distance on R™. The 
admissible values of the constant 7 will be discussed in detail at various instances in the text. 
We impose the initial condition 

u(0, x) = u (x) on {0} x il, (2.8) 

for Uo £ [L 2 (£l)] n . On [0,T] x dil, we impose mixed Dirichlet and Neumann boundary 
conditions as follows. Let i be an index running over the set {1, . . . , n}. For each i the 
boundary d£l is split into dfl = Tf U , with T° being of positive (d — 1) -dimensional 
(Hausdorff) measure. Further, we subdivide dVL = U dQf, where dVLJ := {x € dVL : 
(£?in)(x) < 0} and dflf = d£l\d£l~ are the inflow and outflow parts of the boundary cMl 
for the i-th equation. Finally, we assign 



= 9? 



on if, 

aiVui-n = gf onT^ndnf, (2.9) 
(a, Vm - Bjui) -n = gf on if n 30r , 

for Dirichlet and Neumann data gf e i/ 1//2 (rP), G L 2 (lf), respectively; here and in 
what follows n denotes the unit outward normal vector to dfl. We denote by \i '■ dCt^ K 
the characteristic function of dil^ ■ Upon defining X - := diag(x^, . . . , Xn)> an< ^ : = 
/ — X - , we can write (2.9 1 collectively as 

u = g D on T D and (iVu - X~UB)n = g N on T N , (2.10) 

where g D := {gf , . . . , g°) T and g N := (gf , . . . , 5 ^) T . 

The model problem is completed by interface conditions which are imposed across Tj 
as follows. In view of (|2.5|l, we define the friction coefficients and weights r i; v i ' 2 : Tj — > 



[0, 1] and the permeabilities p { : R 2n — > K", i = 1, . . . , n. The permeabilities are general 
functions of the traces of u from both sides of the interface and thus cover the membrane 
model described in Section \2.\\ . The interface conditions read 

(diVwi - Ui Bj) • n| n i = p^u 1 ,!! 2 ) • (u 2 - u 1 ) - r t {v}u\ + v?uj){Bin)\ n i, on Tj, 
(aiVui - Ui Bj) • n| 2 = p^u 1 ^ 2 ) ■ (u 1 - u 2 ) - r t {v}u\ + v?u?)(Bin)\ Q 2, on Tj, 

where u? := u| nj . nr;j , j = 1,2. Introducing, T 3 = diag(v{, . . . ,v{), j = 1,2, R = 
diag(ri, . . . , r„), and P(u) = (p 1 (u 1 , u 2 ), . . . , p„(u 1 , u 2 )) T , the interface conditions can 
be written in vector notation as 

(AVu - UB) n| n i = P(u)(u 2 - u 1 ) - R(T 1 C/ 1 + T 2 (7 2 )(Bn)| f21 , on IV, 
(AVu-UB)n\ n 2 =P(u)(u 1 -u 2 )-R(T 1 [/ 1 +T 2 ;7 2 )(Bn)| f22 , only (2 ' U) 



We ma ke th e following assumptions on the weights and permeabilities in accordance 

~~ " 12 

For every i = 1, . . . , n, the weights v, ' satisfy, for any x £ Tj, 



with Section 



2.1 



1, ^ 2, ^ 1 J if(An| am )(x)>0, 
"iW+Wi(x) = l, <^ t 2 . (2.12) 

(x) < (x) otherwise. 
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We let p : K — > R n denote the function describing the diffusive flux across the inter- 
face Tj, that is 

p(x\x 2 ) :=P(x 1 ,x 2 )(x 1 -x 2 ) Vx 1 ,x 2 eM", (2.13) 

and assume that p € C 1 ' 1 (M. 2n ) and that its Jacobian p' is uniformly bounded. 

Throughout this work, we shall assume that the above system has a unique solution that 
remains bounded up to, and including, the final time T. 

3. Space discretisation by the discontinuous Galerkin method. 

3.1. Finite element spaces. Let T be a shape-regular and locally quasi-uniform sub- 
division of £1 into disjoint open elements k E 7, such that Tj C U Ke <jdn =: T, the 
skeleton. Further we decompose T into three disjoint subsets T = dft U T mt U Tj, where 
Tin, := r\(c?n U Tj). We assume that the subdivision T is constructed via mappings F K , 
where F K : k — » k are smooth maps with non-singular Jacobian, and k is the reference d- 
dimensional simplex or the reference <i-dimensional (hyper)cube. It is assumed that the union 
of the closures of the elements n E 7 forms a covering of the closure of f2; i.e., Cl = U k6 tK- 

For meNwe denote by V m (k) the set of polynomials of total degree at most m if k is 
the reference simplex, and the set of all tensor-product polynomials on k of degree k in each 
variable, if k is the reference hypercube. Let m K € N be given for each n E 7. We consider 
the /ip-discontinuous finite element space 

V h := {v E L 2 (n) : v\ K o F K E F mK {k), k E 7}, (3.1) 

and setV h := [V,,]". 

Next, we introduce relevant trace operators. Let k + , k~ be two elements sharing an edge 
e := dn + D dn~ C T mt U Tj. Denote the outward normal unit vectors on e of dn + and 3k~ 
by n + and n~, respectively. For functions q : O — > M. n and Q : O — > R nxd that may be 
discontinuous across T, we define the following quantities: for q + := q| K +, q~ := q| K - and 
Q + := Q| K +, Q _ := Q| K - on the restriction to e, we set 

{q}:=l(q++q-), {Q} := ^(Q + + Q"), 

and 

[q] := q + ® n+ + q~ (g) n", [Q] := Q+n+ + Q^n", 

where <Ei denotes the standard tensor product operator, with q ® w = qw T . If e E dn n Tq, 
these definitions are modified as follows: {q} := q + , {Q} := Q + and [q] := q + <S> 
n, [Q] := Q+n. 

Further, we introduce the mesh quantities h:£l— ► R, m:f2— > M. with h(x) = diam k, 
m(x) = m K , if x E n, and the averaged values h(x) = {h}, m(x) = {m}, if x E F. Finally, 
we define ft, max := max ie!) h and h mm := mm xe n h. 

We shall assume the existence of a constant Ca > 1 independent of T such that, on any 
face that is not contained in Tj, given the two elements k, k' sharing that face, the diffusion 
matrix A satisfies 

C-/ < WAW^WA- 1 ^, <C a . (3.2) 

This assumption can be removed using the ideas from [24|, but we refrain from doing so here 
for simplicity of the presentation. 
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The next result is a modification of the classical trace estimate for functions in H 1 (fi 1 U 
Q 2 ) + Vh', see [ 10] for similar results. 

LEMMA 3.1. Assume that tlie mesh T is both shape-regular and locally quasi-uniform. 
Then for v € U fl 2 ) + Vh, the following trace estimate holds: 

2 

X>Mlr 3 <c 1 6(^||V W ||^ + ||h- 1 / 2 M||L)+C2 £ - 1 |kl| 2 , (3.3) 
3=1 /tea" 

for any e > h max and for some constants C\ > and c 2 > 0, depending only on the shape- 
regularity of the mesh and on the domain fl 

Proof. We use the decomposition of v g U £1 2 ) + Vh into a conforming part v c 

and a non-conforming part Vd := v — v c £ Vh- This decomposition is described in ll34l[35l 
for functions in V h \ the extension to i/ 1 (Sl 1 U Q 2 )+V h follows by taking v c £ H 1 (fi 1 Ufi 2 ). 
Using Theorem 2.1(iii) from (35], there exists awje U il 2 ), such that 

E " v c)l!'nn> < CHh 1 ^-!"! [«]||^ nn4) (3.4) 



where 7J Q is the differentiation operator for a multi-index a with \a\ = 0, 1. Hence, (3.4 1 
implies 

for (v c )\w := j = 1,2. 

The triangle inequality implies 

2 2 

ElH^Hr, < 2 E(W^Hr 3 + ll(«-^)|^|iy. (3.6) 

3=1 3=1 



To bound the first term on the right-hand side of ( 3.6 1, we note that v c € H (0 U f2 ), giving 

2 2 

< CE + IKIIn,||V Uc ||ni) 1/2 < C(\\v c \\ 2 + |K||||V*; C ||) 1/2 

3=1 3=1 

<C(e- l \\v c \\ 2 + e\\Vv c f) 1/2 , 

(3.7) 

for any e > sufficiently small. To further bound the right-hand side of (3.7 1, we use the 
triangle inequality for each term, viz., 

ll«c||<||»-«c||+IMI and ||V« C || 2 <2E(I|V(«-« C )|| 2 + ||V«|| 2 ), (3.8) 



in conjunction with ( 3.5 1, to arrive at 

2 

ElM^Hr, <c(e- 1 ||h 1 / 2 H|| 2 mi + e- 1 || U || 2 + e||h- 1 / 2 [i;]|| 2 im +eEl|Vt'|| 2 ) lA> 

3=1 k£7 

< CU-'M 2 + e\\h-^ 2 [v]\\ 2 iM +eE II Vi/"" ' " 



(3.9) 
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using the assumption that e > h max . 

To bound the second term on the right-hand side of (3.6 1, we use the trace estimate on 
each element in conjunction with (3.5 1, viz., 

2 

Ell^-^MIr, <C £ (\\h- l / 2 (v-v c )f K + \\^V(v-v c )f K ) <c\\[v]\\l M . 
j=i Kea:Bnr 3 ^0 

(3.10) 

Noting the assumption e > h max , the use of ( 3.9 1 and ( 3. 10 1 in ( 3.6 1 concludes the proof. □ 

3.2. Space discretization. The discretization of the space variables will be based on a 
discontinuous Galerkin method of interior penalty type for the diffusion part and of upwind 
type for the advection. Special care has to be given to the incorporation of the interface 
conditions. 

More specifically, we shall introduce a interior penalty discontinuous Galerkin (IPDG, 
for short) discretisation of the advection-diffusion operator 

- V- (AVw- WB), (3.11) 

where W := diag(wi, w 2 , ■ ■ ■ , w n ), for w := w 2 , ■ ■ ■ , w n ) T . 
To this end, we define B(\Vh, v/j) to be the IPDG bilinear form 

I (AVw h - W h B) :V Vfl + / ({W h B} + S 3 [w h ]) : [v & ] 

- ^ ({AVw h - W h B} : [v h ] + {AVv h } : [w/J - (S + B)[w h ] : [v fc ]) 

- y ((AVw,, - X+W h B) : (v ft ® n) + (AVv h ) : [w h ® n) - Ew,, • v h ) 
+ / (X+W fe B) : (vfcgn), 



(3.12) 



and define 

JV(w fc ,v h ) := / (P(w fe )[w h ]-(I-R)({T^ h B} + S a [w ft ])):[v /l ]. (3.13) 

Here, S := C cr Am 2 h _1 denotes the discontinuity-penalization parameter matrix with CV > 1 
constant. Furthermore, 23 := | diag(|i?i • n|, . . . , |S„ • n|) and 23j := (T 1 — |l)Bn|ni = 
(T 2 — |l)Bn|na is diagonal with non negative entries. 

REMARK 3.2. We comment on the interface terms. The diffusion term, appearing in N, 
is simply given by J r P(w)[w] : [v^J. This resembles the typical jump stabilisation term 
with the permeability coefficient replacing the discontinuity-penalization parameter, render- 
ing its implementation within a dG computer code straightforward. 

To ensure the coercivity of B, the advective interface term has been split as R = I — (I — 
R), resulting into contributions in both B and N. Indeed, the advective interface contribution 
in B can be recast using the weighted mean {WhB} v := T Wft,B|Qi + T 2 WhB\^2, so that 

{W h BY : [v h ] = ({W h B} + S 3 [w fc ]) : [v h ], (3.14) 

thereby resembling the typical dG upwinding for linear advection problems and, hence, en- 
suring the coercivity of B. 

REMARK 3.3. In this setting, Tn can have non-trivial intersection with both dVL~ and 
9f2^~, thereby extending the discontinuous Galerkin method proposed in A32V . 
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4. Elliptic projection error. The a priori error analysis is based on a (non-standard) 
elliptic projection inspired by a construction of Douglas and Dupont for the treatment of 
non-linear boundary conditions in the context of conforming finite element methods [20]. 

DEFINITION 4.1. For each t G [0, T] we define the elliptic projection -w h G Y h to be 
the solution of the problem: find w/j = Wft(t) G V/,, such that 

B(u-w h ,v h ) + \{u-w h ,v h )+N(u,v h )-N(w h ,v h ) = Vv ft e¥ fe , (4.1) 

for some fixed A > 0. 

The constant A > in the definition above will be chosen large enough to ensure the 



uniqueness of the projection w^, cf. Lemma 4.4 below. 



Next, denoting by § s := HP + V h , s G K, we define the dG-norm on S 1 
|||w||| := ( (WVAVwW* + \\\ Vdiag(V • B)w|| 2 ) + ||VE[w]|£ Durh 

_ V2 

\Vviwl\\*+ ||V^[w]||£ : 



(4.2) 



where ||Q|| 2 := / X)"=i IQiC^)! 2 dx, denotes the Frobenius norm whenever Q is a n x d 
tensor. We assume that j4.2\ is a norm. This is satisfied when standard assumptions on the 
solution in conjunction with the boundary conditions hold on each subdomain, e.g., To n dQP 
has positive (d — 1) -dimensional (Hausdorff) measure for j = 1,2. If the interface manifold 
Fj is not characteristic to the advection field, such hypotheses can be further relaxed. 

For the remainder of this work, we shall make the simplifying assumption that B is such 

that: 



Bi ■ V(v h )i £ V h , for i = l,...,n, (4.3) 

for any function v h := ((i>/t)i, . . . , (f^) n ) T <= V^. We note, however, that this appears not 
to be a genuine limitation in the arguments presented below: ideas on how to circumvent 
this assumption have been presented, e.g., in Il32l l5l. for the case of scalar linear advection- 
diffusion problems. 

The next two results show the coercivity and the continuity of the bilinear form B(-,-). 
Their proofs follow straightforward variations of well-known arguments (see, e.g., flU [32]) 
and are, therefore, omitted for brevity. 

Lemma 4.2. For V/, G V^, there exists a positive constant C coer £ K, such that 

B{v h ,v h ) > C coer |||v h ||| 2 . 

LEMMA 4.3. Let II : [L 2 (0)] n — > denote the L 2 '-orthogonal projection onto Vh- 
For any w G H s , s > 3/2 and v/j G we have 

IS^VfcJI^ac^tllltllllBlllVfclH, 

where r) := w — Ilw and 

INH! := |||r,||| 2 + \\E- 1/2 Wl}\\hur» + ll^Wllr + \\V^{v}fr r (4-4) 
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The next result establishes the well-posedness of the problem ( |4.1| i and relevant approx- 
imation properties. 

LEMMA 4.4. Assume that u G W, s > 3/2 for all t € (0,T]. For A > sufficiently 
large and for h max sufficiently small, the variational problem \4.1\ has a unique solution 
€ Whfor each t € (0, T]. Moreover, the following bound holds: 

Cco^Hlplf + Allpll^lhllll^, (4.5) 

and, if also u t £ HP, then 

C coei .|||p f ||| 2 + A||pJ 2 < |||j7 t ||||,A + IINIil.A (4-6) 

where p := u — w/„ r\ := u — ITu, and 

IIMUl.A :=C c \\\r,\rv+7\\\nf, 

with C c := (4C 2 ont + ZC 2 coer )/C coer . 

Proof. Well-posedness of ( |4. 1 J > is established by proving that the associated mapping is 
strongly monotone on V/j. Using the assumption that p g C 0,1 (]R 2 ™) , we get 

|iV(v,z)-iV(w,z)|< / (|p(v)-p(w)| + C s |v-w)|)|[ Z ]| 
< [ (C|v-w|+C7s|v-w|) |[z]| 



T 3 

2 



(4.7) 



<^£(||(v-w)M 2 3 +||zM| 2 3 ), 

where C p is a Lipschitz constant for the function p, C% > is a constant proportional to 
max HSjiilloo r 3 > an d Cm > is a constant depending on both C p and C3. 

i— l,...,n 

This, in conjunction with the coercivity of the bilinear form B provided by Lemma 4.2 
gives 

B(w h - w fc , v/j - w h ) + A(v,j - w ft , v h - Wfe) + N{v h ,v h - Wh) - N(-w h ,v h - w ft ) 



> CcoerlllVft - W h ||| 2 + A||v h - W,^ 2 - Cg ^ ||(V/, - W/,)| n j 



2 

2 



^ Ccoer 11 |||2 . / \ 2C^CiC2 v., ..■> 
> — 5 — — WA.III +(A-"H 3 jllVh-WfeU 



the last inequality owning to the trace estimate \33\ with e = C coer a m i n /(2C2Ci). Hence 
strong monotonicity is ensured as SOOn as A > C1C2/ Ccoer C^min- 

To show j4.5\ and ( |4.6[ ), we split the error u — w/j = 77 + with 77 := u Liu and 
£ := IIu — Wh- Then, setting = £ in (|4.1|i, we deduce 



I) + A(€, = -B(77, 1) - A(r7, - JV(u, £) + N(w h , £). (4.8) 
Using the coercivity and the continuity of B along with the bound (|4.7|i in (|4.8|l gives 



^ lll^lll^^^^il^ll^ < lll^lli^ -^^M^ll^ ^ C 1 1 ^r? I 1 1 ^ 3 -4- 1 1 I 1 1 > - (4.9) 

^coer ■ , 
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Using (3.3 i with e = C coer a m i n / (4C^Ci) on the last term on the right-hand side of (4.9 1, we 
arrive at 

^ coer |||£||| 2 + ^All^ll 2 < C C \\\V\\\1 + M\V\\ 2 + Cic^WM 2 + ||£|| 2 ). 
Choosing A > l6{C^) 2 ciC2/(C coer a ni i n ), we deduce 

\c coer \u\\\ 2 + \M\a 2 < CcMWl + Im\v\\ 2 - (4.io) 



A triangle inequality already gives ( |4.5) , 

In view of obtaining ( |4.6) , we differentiate ( |4.1[ ) with respect to t and then test with 

(4.11) 



- ^ |(p(u) - p(w h ) + (I — R) ({(17 - W h )B} + <S> 3 \p\) ) 
Using the assumption that the Jacobian p' 6 C 0,1 (IR 2 ") and is bounded we obtain 
J r ~(p(u)-p(w h ) + {I-R)({(U-W h )B} + fy\p])) : fo] 

2 

P '(u) - P '(w ft ))[u t ]| + | P '(w h )[ Pt ]| + c s E(Ip*M + IpM)) IKJI 



3=1 

2 2 

< (E^.ullPMh + C£|| ft M|r 9 )) El^lwllr, 

3=1 3=1 
2 



< 



E(^,ullPl^llr J +2CS||77 t |o 3 '||^ + (C| )U + 3C^)||^ 



2 



it\m\\ri 

3=1 ^ 

2 

< C$,u E OWn* Hr, + lklwllr J + ||€ t | n i||^) J 

3=1 

(4.12) 

with u , u > constants depending on and on ||ut||ioo([o ,t]xsi)- Here and in 
what follows, square brackets are used to denote the argument of a linear operator. 

Applying ( |4.12[ ) to the right-hand side of ( |4. 1 1 1 > and using the coercivity and continuity 
of B yields 

J Ccoer\Mt\\\ 2 + 2 X WZtW 2 < ^Mll^llll + M\Vt\\ 2 
4t 4: ^coer 

2 (4.13) 
+ ^,uE(ll^llr 3 + ll^llr 3 + ll^||^). 

3=1 

As before, using the trace estimate ( 3.3 1 on the last term on the right-hand side of ( 4. 13] > 
with e = C coer a min / (4C| u ci), we arrive to 

^ coe ,|||l t ||| 2 + ^U t \\ 2 < 4Clnt c +3Cler \\\v t \\\l + lM\r) t \\ 2 + JiNIII.a, 
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for any A > 16(max{C^ u , C^}) 2 ciC2/(C coer a m i n ). Now ( |4.6| l easily follows by the trian- 
gular inequality. 

We conclude this section with an L 2 -error bound of the elliptic projection (4.1 1. This is 



obtained by an Aubin-Nitsche duality-type argument, inspired by a construction of Douglas 
and Dupont [20 1 for non-linear boundary conditions. 

The interface operator N given in p,13[ ) consists of a nonlinear component driven by the 
function p(w) = P(w) [w] and a linear component which we can characterise by introduc- 
ing the linear operator L[w] := —(I — R){W / B} tJ + Bj[w]. It is convenient to momentarily 
view N as an operator from indicated with a calligraphic font: 

K: § S*, w >->■ (v h-> J (p(w) + L[w]) : [v] 

Thus the derivative N' is a mapping § — > L(S, §*), where L(S, S*) denotes the linear map- 
pings from § to S* . Therefore the integral 

3>(t,v) := / K'(w e (i,-))(v)^, 



where w 9 := 6u + (1 — 0)w/j, belongs to §* for each i € (0,T), v G S. In particular 
y(t,u(t,-) -w h (t,-)) G §* and 

3>(t, u(t, ■) - w fc (i, •))= / ^'(w e (t, .))(u(i, ■) - w h (t, •)) d# 

9 e (W(w e (t,-)))^ 

o 

= X(u(t,-))-X(w h (t,-)), (4.14) 

using El p.89]. We shall frequently abbreviate CP(i, z(t, •)) by O^z below. 

We assume that there is an s G (3/2,2] such that for all a G [L 2 (fi)] n and /3 G 
[H 1 / 2 (Tj)] 2n there exists a solution £ € H s of the linear dual equation: 

5(v,C) + A(v,0 + (?v,0 = (v,a) + (v,/3) r3 Vv G H 1 . (4.15) 

2 

_f/!/2(r 3 )- (4.16) 

j'=i 

LEMMA 4.5. Assume that the hypothesis of Lemma \4.4\ and ( |4.15| > vwf/i ( |4.16[ ) hold true. 
For A > sufficiently large, for /i max sufficiently small, the following error bounds holds: 

IIpII < c(i + ^L x A) 1/2 ^r a x IIMIkAi (4-i7) 

If in addition the Hessian p" is uniformly bounded and u, u t G W 1,oo ([0, T] x Tj) f/zen 

|| Pt || 2 < C(l + < ax A) 1/2 ^ ax (UMII S , A + IIMIM- (4-18) 
The constant C depends only on Ca and the shape-regularity of the mesh. 
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Proof. Let z solve ( |4.15| ) with a = p and (3 = 0. Owing to the adjoint consistency 
of the symmetric interior penality bilinear form B as well as of CP, the dual solution z also 
satisfies ( |4.15| l for v e V/i, thus allowing us to test in ( |4.15| l with v = p = u w/j and get: 

||p|| 2 = B(p,z) + A(p,z) + (CPp,z) 

= B(p, z) + \{p, z) + N(u, z) - N(w h , z) (4. 19) 

= B(p, n z ) + X(p, V z ) + N(u, r] z ) - N(w h ,r, z ), 

with r/ z = z IIz. The last equality follows from ( |4.14| ) and the definition of the elliptic 



projection (4.1 1 



We now bound each term on the right-hand side of ( |4.19| l. Observe that 

|B(p, J 7 z )|<qi|p|||(|||^|||| + E(_ max h7 1/2 ml^\\ri z \\l) 1/2 

ket 1 (4.20) 
< C|||p|||(|h*|||| + A||^|| 2 ) 1/2 <C|||p||||||77 z ||| B , A , 

for A big enough. The nonlinear interface term in ( |4.19| l is bounded as in the proof of 
Lemma |44j yielding 

2 2 

|7V(u,, 7 2 )-iV(w, I ,r ? z )|<4C|(^||p| fP ||^) 1/2 (^||^| 03 ||2 3 ) 1 / 2 

3=1 3=1 

(4 21) 

< (C coer |||p||| 2 + A||p|| 2 ) 1/2 (C coer |||r7 z ||| 2 + A||rf|| 2 ) 1/2 , 



using once again (3.3 1 with e = C coer a min /(4C^c 1 ), for A > 16(C^) 2 CiC 2 / (C coer a min ) 



Using the above inequalities and applying the bounds on the elliptic projection error 



given in Lemma 4.4 we obtain 

l|p|| 2 < C|IMII a ,Jto'lll Bl A (4-22) 

for A big enough. 

The term in rf can be bounded using standard approximation estimates for the error of 
the (orthogonal) L 2 -projection (see, e.g., (48j) yielding 

^ (S_1) (1 + KX) ||z||| s(k) . (4.23) 



Further, inserting (4.23 i into (4.22 1 and then using ( |4.16| ) we get 



2 

||p|| 2 < C(l + hl^XY^h*-^ |||r,||| S)A J2 Mn-m 

3 = 1 

<a(i + ftLxA) 1/a /cilhlll s ,xllp|l. 

thus yielding ( |4.17| >. 



(4.24) 
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Testing in ( |4.15| ) with v = p t gives, for each t, 

\\ Pt (t)\\ 2 = B{p t ,i) + A( Pi! z> + (3>p t> 2> 

= B(p t ,z) + X( Pt , z) + (dt?(t, p(t)),z) - (d t 9, z) 

= Bfa,*) + A(p t ,i) + (d t (N(w, •) - N(w h , .)),*)-(%?, z) 

= B(p t , V s ) + A(p t , rf) + (d t (JV(w, •) - AT( Wft , •)) , r, 2 )-(d t $, z). 



(4.25) 



with r/ z = z — IIz. The last equality follows from differentiating with respect to time the 
definition of the elliptic projection (4.1 1. 

The first three terms on the right-hand side of ( |4.25[ ) can be bounded by the same ar- 
gument used above. In particular, following the argument in ( |4.12| i and then applying the 
bounds on the elliptic projection error given in Lemma |4~4] yields 



\dt (N(u, V z )~N(^ h ,rj B ))\ 



<4^,u(£(IIPM& + IIPtlrfJ) 7 (£> If-IU- 
3=1 j'=i 

<C(\\\r,\\L x - 



2 \l/2 



(4.26) 



As for the last term, we proceed as follows. Spelling out the definition of 7, we have 

(%, P (t)),z) = 



^(w e )(p)^, z ; 

^((pCw^+L)'^],!!])^ 

(p'(w e )[p])+L[p]):[z]W0. 



Thus, using the assumption that the Hessian p" i s bo unded, the embedding of iP(fP) into 
L 00 ^), j — 1, 2, and then using ( |4.16| l, Lemma 3.1 and ( |4.6| l, we get: 



l(a t y| (tiP(t)) ,z)| 



< 







(p"(w e )[a t w e ,p]):[z]j^ 
|p"(w e )P t w e ||p||[z]|) dfl 



(4.27) 



<<?(£ Noo.ro) / (|9 t p| + |9 t w A |)|p| 

<C u ||p t ||||p|| H1 / 2( r 3 ).. 

Hence we are left in need of an estimate of ||p|| i ji/2(r/ 3 )«- This can be obtained by the 
following duality argument. 

Let z be the solution ( |4.15[ ) with a = and (3 = 5(p). Here S is the duality map [51, 
HB, p. 860] between tf 1/2 (T;,) and H 1 / 2 ^^* such that \\8(p)\\ H i/2 (r:i) = \\p\\ H ^^r 3 )' 
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and (p,5(p))r 3 = ||p||^i/2( rj )*- Testing in ( |4.15| > with v = p yields 

||p|lHV=(r 3 )* = (P' 5 (P)>r 3 

= J B(p,z) + A(p,z> + (Tp,z), 

= B(p, z) + \{p, z) + N(u, z) - AT(w ft! z) 

= B(p, rf ) + A<p, tj*) + JV(u, rf ) - AT( W/l , ), 

with r/ z = z ilz. 

Now a bound for ||p||#i/2(i\,)* can be derived by following the same steps used above 
to get ( |4.24[ ), yielding: 

2 

||p|llri/2 ( r 3 ). < C(l + /ii ax A) 1/2 ^r a x \\\v\h,x J2 H 4 IU'(no 

3=1 (4.29) 

< C(l + /iL x A) 1/2 ^ ax MhJ^Wm/^), 

having used ( |4.16[ ) once again. Finally, by inserting in ( |4.29| l the definition of S(p) we con- 
clude that 

\\p\\ H U2 [T3) , < C(l + h 2 max \)V 2 h s -^ ||MU S , A . (4.30) 

Using ( |4.30| > in ( |4.27 1 we bound last term on the right-hand side of ( |4.25| l. Now, recall- 
ing \A.26) , the bound ( |4.18| ) easily follows from \A.25\ . 

5. DG method for the parabolic system and its error analysis. The above discussion 
motivates the introduction of the following IPDG-in-space method for the system \2.6) , ( |2.8| >, 
plQl , and (fZlT). 

For t = 0, let u fc (0) = w/,(0). For f G (0, T], find u h = u h {t) G V/, such that 
((u h )t,Vh)+B(u fc ,v h )+^(ufc,Vfc) + <P(uh) s Vh)=l(v/i), for all v h G V h , (5.1) 
where 

J(vj,) :=- / ((g D ®n) : (4Vv A ) + (rG D B) : (v fe (g>n) -£g D - V/l ) + / g N -v h , (5.2) 

noting that (q ® n) : (v ® n) = q • v, for q, v G R™, and having denoted Go := 
diag(5D,...,5g). 

We shall make use of the following result to treat the nonlinear reaction term. 

LEMMA 5.1. If d = 2, let q G [1, oo) anrf f/d = 3, let q G [1, 6]. 77zen, f/zere erafs a 
constant Cpf > 0, depending only on the geometry of the subdomains , j G {1, 2} one/ o/ 
?/ie Dirichlet boundary, such that, for all v G S , we have 

\\v\\l<C v ?m^{l,a^}\\\v\f. (5.3) 

Moreover, for < "/ < 2, if d = 2, and for < 7 < 4/3, i/d = 3, we /lave 

||«||^ < C PF max{l, am i n }||t;|n||«||| 2 . (5.4) 



Finally, for d = 2 and 7 = 2 f which corresponds to q = 00), ( 5. 5 1 ant/ ( 5.4 1 hold for v G V^, 
w/f/z Cpf := C| log(mino h) | /or some constant C > depending only on Q l . 



16 



A. CANGIANI, E. H. GEORGOULIS AND M. JENSEN 



Proof. From the assumptions on the topology of the subdomains, we can apply Theorem 
3.7 from [39], to deduce 

|M|^ <Cmax{La m i n }||M|| 2 , 

with the boundary contribution taken as dfl^Tj, for l<<7<ooif<i = 2 and for 1 < q < 6 
if d = 3. For the case q — oo if d = 2, we make use of the standard inverse estimate 
IMIL, K < dogl^HMIff^K)- Setting p = 2/ 7 , for 1/p + l/q = 1, Holder's inequality 
implies ( |5.4| >. I_l 

We are now ready to prove a bound for the difference between the elliptic projection w; t 
and the dG approximation u/j. 

THEOREM 5.2. Consider the notation of Lemma \4.4\ and let 7 as in Lemma \5J\ Assume 
thatu e L 2 (0,T;HP) n L°°(0,T x fi), s > 3/2, u t £ L 2 (0, T; L 2 (Vt)), and that ft, max is 
small enough. Then, we have 

(0,T;L 2 (a)) + Ccoer 1 1#| I L 2 (0 : T;S) - 4(5 2 e CT , 

where 6 := Wh — U/j, C 

s 2 (t) 



C(||u||oo,[o,T]xn,A), and 



C'llp|| 2 + C(||p||^ + A- I ||pj 2 )di. 



(5.5) 



Proof. Without loss of generality assume that 7 > 0; the case 7 = follows by a simple 
modification of the argument presented below. 

Let e := p + 8, with p := u — w/, and 9 := w/j — u^. We note that the continuity of w/j 
in the time variable is implied by the well-posedness of the elliptic projection problem, and 
that of U/j by a standard local existence argument near on the resulting system of ordinary 
differential equations. Hence, 6 is continuous in [0, T] with 6(0) = 0. 

Orthogonality implies 

(d t e, 9) + B(e, 9) + N(u, 9) - N(u h , 9) + (F(u) - F(u h ),9) = 0. 
From ( |4.1| ), we deduce 

^||0|| 2 + B(9, 9) = (F(u h ) - F(u), 9) +N(u h , 9) - N(w h , 9) + (Xp - p t ,9). (5.6) 
Making use of the inequality 

for a, b > 0, the non-linear reaction term can be bounded as follows: 
|<F(u)-F(u„),0)| <C f (l + \up + \u h p)\u-u h \\9\ 

<c /(i+|uniu-u fc ||0|+c / iu- U7l r+ i |0| ( 5 - 7) 

Jn Jn 
<C(u)(||p|| 2 + ||0|| 2 )+C(||p||^ + ||0||^ 2 ). 
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Also, using the regularity of p and (3.3 I, we have 



\N(u h ,e)-N^ h ,e)\<C^Y l W y\\h < ^C cocr |||0||| 2 + ^||<9|| 2 , (5.8) 



choosing e and A as in the proof of Lemma 4.4 
We also have 



|(Ap-p t ,0)|<^||p|| 2 + ^|| Pt || 2 + A||0|| 2 . (5.9) 



Lemma |5.1| implies 



Wew^tl < c PF m^{i,a-l}\\er\\\e\\\ 2 . (5.10) 

Using these bounds in (|5.6|l, along with the coercivity and continuity bounds from Lemmas 



4.2 and 4.3 integrating the resulting inequality with respect to t (and multiplying by 2) be- 



tween and t, we arrive at 

||0(r)|| 2 + C coer f |||0||| 2 
Jo 

<s 2 + c [ T \\e\\ 2 + Ce S ssu P \\8(tw [ T \\\e\\\ 2 (5 . n) 

Jo te[o,r] Jo 

<5 2 + cf \\0\\ 2 + c(esssup\\9(t)\\ 2 + C coer f |||0|] 
Jo ^ te\o,T] Jo 



1+7/2 



for S 2 and C as in the statement of the theorem, and C := C max{l, (C coer a m i n ) 1 }. From 
the assumption that the mesh-size h niax is small enough, we can have 

which implies (4(7<5 2 e CT ) 1+7 / 2 < 5 2 . We now consider the set 

/:= {re [0,T] :csssup||6»(<)|| 2 + C coer f |||0||| 2 < 4<5 2 e (5T }, 
te[o,r] Jo 

which is non-empty and closed due to the continuity of with respect to the time variable, as 
8(0) = 0. We set r* = max / and we suppose that r* < T. Hence, for r < r* , we have 

esssu P ||0(t)|| 2 + C coer f|||0||| 2 <2<5 2 + <7 f \\6\\ 2 . 
te[o,r] Jo Jo 

Gronwall's Lemma then implies 

\\0(r*)\\ 2 + C coer f \\\6\\\ 2 < 25 2 e CT , (5.12) 
Jo 

setting t = t* , which contradicts the hypothesis t* < T, due to the continuity of the left- 



hand side of (5.12 1. Hence, / = [0, T] and the result holds. I_l 



COROLLARY 5.3. Let < 7 < 2 d = 2, and < 7 < 4/3 // d = 3. With the 
assumptions of Lemma 4.4 and of Theorem 5.2 and u| K € H l (§,T\ [_ff" +1 («;)]"), k K > 1, 



k6X we /lave 

||u - Uh||i«»(o,Tii=(n)) + c 'co e r||u-Uh||| a(0)T . s) < C£((0, T], h, u, V h ), (5.13) 
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for C independent ofh and 

£((0, T],h, u,W h ) := £ [ T hl°« (|u|f Hfcre+1(fi)]n + \vL t \\ HkK+l{K)]n ), (5.14) 



for s K = min{m K , k K }. Moreover, with the additional assumptions of Lemma 4.5 we have 



llu - Uh||£co ( o ir;L9(n) ) < Ch s -^£((0,T},h,u,Y h ), (5.15) 
for s G (1,2] dictated by the reg ular ity of the exact solution of the dual problem ( |4.15| l. 



Proof. We begin by using (5.3 1 to deduce the estimates ||p|| 7 + 2 — C|IIpIII 7+2 > ||p| 



< 

C|||p|||, and ||p t || < C|||p t |||. An application of Lemma 4.4 on the resulting norms, along 
with standard approximation estimates for the error of the (orthogonal) i 2 -projection (see, 
e.g., 11321 ) already implies ( |5.13| l. 

For ( |5.15| l, we set e = min{2, kj/(s — 1)}. Then, we have the bound 

\\ P \\;t 2 2 < iipiMipii^/e+D < ciipii 2 - £ nipiir +e , 

which implies (for this choice of e, since then 7/e + 1 < 2) that 

\\p\^Xl<Ch s -^e.((o,T],h,u,w h ). 

Lemma 4.5 and standard approximation estimates now imply ( |5.15| l. □ 



REMARK 5.4. We note that the bounds in Theorem \5.2\ and, correspondingly, in Corol- 
lary \5.3\ do not require any global mesh quasi-uniformity assumptions. This result also holds 
on domains without internal interfaces as in the setting of H40V . There, a different continu- 
ation argument is used for deriving a priori bounds which requires the mesh to be globally 
quasiuniform (for non-symmetric spatial operators), albeit for a larger range 0/7 than the 
one considered in the present work. 

REMARK 5.5. It is interesting to note that, due to the careful use of L 2 -projection 
operators in conjunction with assumption ( |4.3| l, the constant in \5.\3\ depends on the Peclet 
number only through the control of the nonlinear interface conditions. Thus, the present 
bound produces error control that is Peclet number independent for the dG method applied 
to a single domain problem. 

It is possible to show optimal error estimates with less restrictive assumptions than ( |2.7| i 
on the growth of the reaction term. Indeed, assuming only that F is locally Lipschitz, an op- 
timal a priori bound can be proven, subject to certain conditions on the mesh. This argument 
is motivated by ideas presented in [2, 25 1 for different problems. 

To this end, consider F l '■ R™ — > R n satisfying 

|F L (x)-F L (y)| <C L |x-y|, (5.16) 

such that F(x) = F L (x), for all x G 1" with |x| < L := 2max < t < T IKi)^. This 
implies, in particular, that Fl(u) = F(u). 

THEOREM 5.6. Consider the notation of Lemma\4.4\and the assumptions of Lemma\4.5\ 



and of Theorem 5.2 Assume also that u| K € H 1 (0,T;[H k ' t+1 (K)] n ), k K > 1, K € land 



that the mesh T is fine enough so that 

hJ n K~lZ((0,T},h,u,W h ) (5.17) 

is small enough ( cf. ( |5.14[ ) for the definition of£. ). Assume, finally, that F is locally Lispschitz. 
Then, we have 

\\vL-u h \\ L oo i0iT . LH a)) < Ch s -^E((0,T],h,n,Y h ), (5.18) 
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with C independent of In. 

Proof. Assume initially that the locally Lipschitz continuous F is replaced with the 
globally Lipschitz continuous Fx from ( |5.16| l, and consider the modified initial/boundary 
value problem described in Section 2.3 with F replaced by F^. Noting that Fl(u) = F(u), 
we conclude that the analytical solution of the modified and of the original problem coincide. 
Let um denote the numerical solution of the modified problem by the dG method \5.\\ with 
F replaced by F L . 

Let = p + 6l, with p = u and 0l ■= w/ t — u^. Orthogonality implies: 

(d t e L ,0 L ) + B{e L ,6 L ) + N(u, 6 L ) - N(u Lh , 9 L ) + (F(u) - F L (u Lh ), 9 L ) = 0. 



Owing to (4.1 1, this gives 

^ll^ll 2 + B{6 Ll 6 L ) = (F L (u Lh ) - F(u),0 L ) 

+ N(u Lh , 6 L ) - N{w hl 8 L ) + (Ap - Pt ,9 L ). 

The terms involving the semilinear form N(-,-) can be bounded in a completely analogous 
fashion to (5.8 1, while the last term on the right-hand side of (5.19 1 can be treated as in (5.9 1. 
Since F l (u) = F(u), the reaction term can be bounded as follows: 

|(F(u)-F L (u Lh ),0 L )| <C L J \u-u Lh \\0 L \ < *C L (\\p\\ 2 +3\\e L f). (5.20) 

Hence, ( |5.19[ ) gives 

\\0 l (t)\\ 2 + C coer ( \\\e L \\\ 2 <S 2 +3(C L + A) f \\8lW 2 , (5.21) 
Jo Jo 

with 8l(t) = Jo T (C L + A)||p|| 2 + A- 1 ||p t || 2 di, noting that u Lh (0) = u h (0). Gronwall's 
Lemma then implies 

||0 L (T)|| 2 + C coer f |||^| 2 <<S|e 3 ^ +A > T (5.22) 
Jo 



Hence the triangle inequality implies 

||ei|Uoo (0 ,r;i3(fi)) < 8 L e z{CL+ ^ T l 2 + ||p|| L ~(o,T;L>(n)) 
<Ch s -^E((0,T),h,u,V h ), 



(5.23) 



using Lemma 4.5 



We shall show that under the mesh assumption ( |5.17[ ), the bound ( |5.23| l also holds for 
u — Uft . To this end, consider the standard nodal interpolation operator 3 K : H s (k) DC(k) — > 
P mic , (see, e.g., IfTBI for the scalar version), satisfying 

\v - 3 K v\ HS(K) < Ch s - j \v\ H s (K) , (5.24) 

for < j < s and s > 2, and 

||« - 3 K w||oo, K < Ch*\v\ w a,oa( K ). (5.25) 
Let also := J re i>i, for v = (vi, . . . , v n ) £ [H s (n) D C{R)\ n . Then we have 

max IIul^IIoo < max - JuHoo + max ||u — "Ju\\oo + max Hu^, (5.26) 

0<t<T 0<t<T 0<t<T 
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The second and third terms on the right-hand side of ( |5.26| > can be bounded using (5.25 1 and 
the definition of L, respectively, giving 



0<t< X T I'" 1 '' 1 " 00 " (" UL ' 1 ~ 3U ^°° + C {^ fc «l U lw=-~(«)) 1 ) + \- ( 5 ' 27) 

For the first term on the right-hand side of ( j5.27| i, a standard inverse estimate implies 

||u ifc - JulU < C'Y, \\ulh - M\* <C<Y, h -~ 2 (ll e ilU + ll u - M\k)- 
Therefore, in view of ( |5.23| l and ( |5.24| >, we deduce from ( |5.27| i the bound 



max Hu^lU <C/ lm |^max£((0,T],h,u,V / + c(^^ |u| 2 i/2 =o(K ^ _ 

kGT 



1/2 L 

2"' 

(5.28) 

Choosing h such that the first two terms on the right-hand side of ( |5.28| l are bounded strictly 
by L/2, one finds ulh = u^, thereby concluding the proof. □ 

6. Numerical examples. In view of the numerical tests, we introduce a fully discrete 
discretisation for the system (2.6 1, (2.8 1, (|2.1 The time discretization is based on a 2nd- 



order linearly implicit method analysed in ||2] [TJ . Let us write the semidiscrete DG formula- 
tion ( |5.1| l in matrix notation as 

MU t = LU + F(U) (6.1) 

where L collects all linear terms and F the nonlinear terms in ( |5.1| l. Treating the linear terms 
with the second-order Adams-Moulton method (trapezium rule) and the non-linear terms with 
the second-order Adams-Bashforth method yields the following fully discrete method (AB2- 
AM2): 

MU n+1 = MU n + k(6AU n+1 + (1 - 8)AU n ) + ^(3F(C/ n ) - F{U n ^)). 

Convergence test. We test the validity of the IPDG method and above error bounds on 
a system of two equations for which the exact solution is known. 

Let the domain = [— 1, l] 2 be subdivided into two subdomains interfacing at x = 0; 
thus n 1 = [-1,0] x [-1,1] and O 2 = [0,1] x [-1,1]. We set T D = {±1} x [-1,1] and 
Tn = (—1, 1) x {±1} and impose homogeneous Dirichlet and Neumann boundary conditions 
on Td and Tn, respectively, as shown in Figure 6.1 (left). For t > we consider the system 
of two advection-diffusion equations 



u t — Aw — u x = f u 



v(l-v) inO 1 



-v 



in ft 2 , (6.2) 



v t - Au - v x = f v + u in n 1 U O 2 , 



and accordingly construct the forcing terms f v ' v : [0, 1] x ft —> M in order to yield as exact 
solution 



it 



/ cosi\ i y 2 -i) 2 j 4a; (1 4- x) infi 1 , 
v J ~ \ sinf ) C X | (-4x 3 + 3a; + 1) inO 2 . 
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4t 




u = 



FIG. 6.1. Convergence test. Solution domain and boundary conditions (left). First component of the exact 
solution at the final time t = 1 (right). 



The first component of the solution at time t = 1 is shown in Figure 6.1 (right). The func- 
tions in ( |6.3[ > are compatible with the interface conditions ( |2.11| i with respect to the interface 
parameters 

P = diag(3,3), Y 1 = diag(l, 1), T 2 = diag(0, 0), R = diag(l,l), 

which are therefore used to close problem ( |6.2[ ) with conditions ( |2.1 1| >. We tested the con- 
vergence rate of our method under space discretisation refinement. A fixed time step of size 
.5 x 10~ 3 is used throughout while an initial uniform square 4x4 mesh is uniformly refined. 
The value C a = 10 was used throughout for the discontinuity-penalization parameter costant. 
Numerical results are reported in Table |6]in the cases of bilinear and biquadratic spatial dis- 
cretizations. The predicted error rates of convergenece are confirmed in both the L 2 (0, 1; S) 
andi 00 ^, l;L 2 (fi)) norms. 

Advection-dominated test. With this numerical example we test the robustness of the 
method in the advection-dominated regime. We consider once again the domain Q, = [—1, l] 2 
subdivided into the subdomains ft 1 and f2 2 already considered in the previous test. In fi, we 
solve a scalar equation with no reaction and the constant diffusion a\ = 10~ 2 and the advec- 
tion field B\ = (0.5,0.5). On the boundary dft, we set homogeneous Neumann boundary 
conditions. The parameters 

P = 2/10, ^ = 5/6, T 2 = l/6, R = 6/10, 

are used in the interface conditions ( |2.1 In this case, transfer across the interface is mainly 
advection-driven. Further, setting the friction coefficient to less than one models the case in 
which the interface acts as a filtering wall on the advected quantity; hence a boundary layer 
in the upwind subdomain in the proximity of the interface is expected. 

We solve the problem on a uniform 16 x 16 mesh using bilinear elements. Such mesh is 
not fine enough to resolve the layer forming in the proximity of the interface where the solu- 
tion is also discontinuous, see Figure [5~2| Nevertheless, the numerical solution is stable, and 
the expected behaviour of the solution is accurately captured, as we can see by comparison 



with the solution obtained with a layer resolving 64 x 64 mesh and shown in Figure 6.3 



7. Concluding remarks. A dG method for the numerical solution of nonlinear inter- 
face problems modelling mass transfer through semi-permeable membranes is presented and 
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# cells 


#dofs 


L a (0,l;S) 


L°°(0,l;i 2 (r!)) 


m= 1 


4 


32 


1.579e+01 




2.498e+00 




16 


128 


7.617e+00 


1.05 


8.293e-01 


1.59 


64 


512 


3.615e+00 


1.08 


2.400e-01 


1.79 


256 


2048 


1.736e+00 


1.06 


6.497e-02 


1.89 


1024 


8192 


8.475e-01 


1.03 


1.693e-02 


1.94 


4096 


32768 


4.182e-01 


1.02 


4.324e-03 


1.97 


m = 2 


4 


72 


4.781e+00 




3.812e-01 




16 


288 


1.013e+00 


2.24 


5.422e-02 


2.81 


64 


1152 


2.282e-01 


2.15 


7.330e-03 


2.89 


256 


4608 


5.480e-02 


2.06 


1.240e-03 


2.56 


1024 


18432 


1.349e-02 


2.02 


1.743e-04 


2.83 


4096 


73728 


3.427e-03 


1.98 


2.259e-05 


2.95 


Table 6.1 



Convergence test. Errors and convergence rates under uniform mesh refinement. A fixed time step of size 
.5 X 10 - 3 was used. Bilinear (above) and biquadratic (below) dG in space discretizations. 




FIG. 6.2. Advection-dominated test. Snapshots of the solution computed on a uniform 16 X 16 mesh using 
bilinear elements: the initial condition ( top-left) followed by the solution at time intervals of 0.5. 



a priori error bounds are shown under typical regularity assumptions. The good performance 
of the method is highlighted through numerical experiments. A number of extensions of 
the presented results can be made with modest modifications. For instance, /ip-version error 
bounds can be shown and more general convection coefficients B can be treated. We re- 
frained from doing so in the interest of simplicity of the presentation. Interesting directions 
of further research are the consideration of the variational crimes due to inexact representation 
of the interface manifold, using, e.g., unfitted finite elements ll7l [T5ll or the related approach 
in flTt and the treatment of more general interface nonlinearities. These will be considered 
elsewhere. 
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FIG. 6.3. Advection-dominated test. Snapshots of the solution computed on a uniform 64 X 64 mesh using 
bilinear elements: solution at t = 0.5 (left) and t = 1 f right) corresponding to the second and third plots in the first 
row of Figure\6.2\ respectively. 



REFERENCES 



[1] G. AKRIVIS AND M. CROUZEIX, Linearly implicit methods for nonlinear parabolic equations, Math. Comp., 

73 (2004), pp. 613-635 (electronic). 
[2] G. AKRIVIS, M. CROUZEIX, AND C. MAKRIDAKIS, Implicit-explicit multistep methods for quasilinear 

parabolic equations, Numer. Math., 82 (1999), pp. 521-541. 
[3] G. D. AKRIVIS AND V. A. DOUGALIS, Finite difference discretizations of some initial and boundary value 

problems with interface, Math. Comp., 56 (1991), pp. 505-522. 
[4] D. N. ARNOLD, An interior penalty finite element method with discontinuous elements, SIAM Journal on 

Numerical Analysis, 19 (1982), pp. 742-760. 
[5] B. AYUSO AND L. D. MARINI, Discontinuous Galerkin methods for advection-diffusion-reaction problems, 

SIAM J. Numer. Anal., 47 (2009), pp. 1391-1420. 
[6] I. BABUSKA, The finite element method for elliptic equations with discontinuous coefficients, Computing 

(Arch. Elektron. Rechnen), 5 (1970), pp. 207-213. 
[7] J. W. BARRETT AND C. M. ELLIOTT, Fitted and unfitted finite-element methods for elliptic equations with 

smooth interfaces, IMA J. Numer. Anal., 7 (1987), pp. 283-300. 
[8] S. BARTELS, A posteriori error analysis for time-dependent Ginzhurg-Landau type equations, Numer. Math., 

99 (2005), pp. 557-583. 

[9] C. BERNARDI, Y. MADAY, AND F. RAPETTI, Basics and some applications of the mortar element method, 

GAMM-Mitt., 28 (2005), pp. 97-123. 
[10] A. BUFFA AND C. ORTNER, Compact embeddings of broken Sobolev spaces and applications, IMA J. Numer. 

Anal., 29 (2009), pp. 827-855. 
[11] E. BURMAN AND P. ZUNINO, A domain decomposition method based on weighted interior penalties for 

advection-diffusion-reaction problems, SIAM J. Numer. Anal., 44 (2006), pp. 1612-1638 (electronic). 
[12] Z. CAI, X. YE, AND S. ZHANG, Discontinuous Galerkin finite element methods for interface problems: a 

priori and a posteriori error estimations, SIAM J. Numer. Anal., 49 (2011), pp. 1761-1787. 
[13] F. CALABRO AND P. ZUNINO, Analysis of parabolic problems on partitioned domains with nonlinear condi- 
tions at the interface. Application to mass transfer through semi-permeable membranes, Math. Models 

Methods Appl. Sci., 16 (2006), pp. 479-501. 
[14] A. CANGIANI AND R. NATALINI, A spatial model of cellular molecular trafficking including active transport 

along microtubules., Journal of Theoretical Biology, 267 (2010), pp. 614-625. 
[15] Z. CHEN AND J. ZOU, Finite element methods and their convergence for elliptic and parabolic interface 

problems, Numer. Math., 79 (1998), pp. 175-202. 
[16] P. G. ClARLET, The finite element method for elliptic problems, vol. 40 of Classics in Applied Mathematics, 

Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. Reprint of the 1978 

original [North-Holland, Amsterdam; MR0520174 (58 #25001)]. 
[17] B. COCKBURN AND C. DAWSON, Approximation of the velocity by coupling discontinuous Galerkin and 

mixed finite element methods for flow problems, Comput. Geosci., 6 (2002), pp. 505-522. Locally con- 
servative numerical methods for flow in porous media. 
[18] D. A. Dl PlETRO AND A. ERN, Mathematical aspects of discontinuous Galerkin methods, vol. 69 of 

Mathematiques & Applications (Berlin) [Mathematics & Applications], Springer, Heidelberg, 2012. 
[19] V. DOLEJSI, M. FEISTAUER, V. KUCERA, AND V. SOBOTIKOVA, An optimal L°° (L 2 )-error estimate for 



24 A. CANGIANI, E. H. GEORGOULIS AND M. JENSEN 

the discontinuous Galerkin approximation of a nonlinear non-stationary convection-diffusion problem, 

IMA J. Numer. Anal., 28 (2008), pp. 496-521. 
[20] J. DOUGLAS, Jr. AND T. DUPONT, Galerkin methods for parabolic equations with nonlinear boundary 

conditions, Numer. Math., 20 (1972/73), pp. 213-237. 
[21] M. DRYJA, On discontinuous Galerkin methods for elliptic problems with discontinuous coefficients, Comput. 

Methods Appl. Math., 3 (2003), pp. 76-85 (electronic). Dedicated to Raytcho Lazarov. 
[22] A. ERN, I. MOZOLEVSKI, AND L. SCHUH, Discontinuous Galerkin approximation of two-phase flows in 

heterogeneous porous media with discontinuous capillary pressures, Comput. Methods Appl. Mech. 

Engrg., 199 (2010), pp. 1491-1501. 
[23] A. ERN AND J. PROFT, Multi-algorithmic methods for coupled hyperbolic-parabolic problems, Int. J. Numer. 

Anal. Model., 3 (2006), pp. 94-114. 
[24] A. ERN, A. F. STEPHANSEN, AND P. ZUNINO, A discontinuous Galerkin method with weighted averages 

for advection-diffusion equations with locally small and anisotropic diffusivity, IMA J. Numer. Anal., 29 

(2009), pp. 235-256. 

[25] X. FENG AND O. A. KARAKASHIAN, Fully discrete dynamic mesh discontinuous Galerkin methods for the 
Cahn-Hilliard equation of phase transition. Math. Comp., 76 (2007), pp. 1093-1117 (electronic). 

[26] M. H. FRIEDMAN, Principles and Models of Biological Transport, Springer, 2008. 2nd ed. 

[27] V. GlRAULT AND B . RIVIERE, DG approximation of coupled Navier-Stokes and Darcy equations by Beaver- 
Joseph-Saffman interface condition, SIAM J. Numer. Anal., 47 (2009), pp. 2052-2089. 

[28] G. GUYOMARC'H, C.-O. LEE, AND K. JEON, A discontinuous Galerkin method for elliptic interface prob- 
lems with application to electroporation, Comm. Numer. Methods Engrg., 25 (2009), pp. 991-1008. 

[29] P. HANSBO, Nitsche's method for interface problems in computational mechanics, GAMM-Mitt, 28 (2005), 
pp. 183-206. 

[30] B. HEINRICH AND S. NlCAISE, The Nitsche mortar finite-element method for transmission problems with 

singularities, IMA J. Numer. Anal., 23 (2003), pp. 331-358. 
[31] J. S. HESTHAVEN AND T. WARBURTON, Nodal discontinuous Galerkin methods, vol. 54 of Texts in Applied 

Mathematics, Springer, New York, 2008. Algorithms, analysis, and applications. 
[32] P. HOUSTON, C. SCHWAB, AND E. SULI, Discontinuous hp-finite element methods for advection-diffusion- 

reaction problems, SIAM J. Numer. Anal., 39 (2002), pp. 2133-2163 (electronic). 
[33] G. KANSCHAT AND B. RIVIERE, A strongly conservative finite element method for the coupling of Stokes 

and Darcy flow, J. Comput. Phys., 229 (2010), pp. 5933-5943. 
[34] O. A. KARAKASHIAN AND F. PASCAL, A posteriori error estimates for a discontinuous Galerkin approxi- 
mation of second-order elliptic problems, SIAM J. Numer. Anal., 41 (2003), pp. 2374-2399 (electronic). 
[35] , Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic problems, 

SIAM J. Numer. Anal., 45 (2007), pp. 641-665 (electronic). 
[36] A. KATCHALSKY AND P. CURRAN, Nonequilibrium Thermodynamics in Biophysics, Harvard University 

Press, 1981. 

[37] O. KEDEM AND A. KATCHALSKY, Thermodynamic analysis of the permeability of biological membrane to 

non-electrolytes, Biochimica et Biophysica Acta, 27 (1958), pp. 229-246. 
[38] D. KESSLER, R. H. NOCHETTO, AND A. SCHMIDT, A posteriori error control for the Allen-Cahn problem: 

circumventing Gronwall's inequality, M2AN Math. Model. Numer. Anal., 38 (2004), pp. 129-142. 
[39] A. LASIS AND E. SULI, Poincari-type inequalities for broken Sobolev spaces, Oxford University Computing 

Laboratory, Tech. report NA-03/10 (2003). 
[40] A. LASIS AND E. SULI, hp-version discontinuous galerkin finite element method for semilinear parabolic 

problems, SIAM Journal on Numerical Analysis, 45 (2007), pp. 1544-1569. 
[41] J. Li, J. M. MELENK, B. WOHLMUTH, AND J. Zou, Optimal a priori estimates for higher order finite 

elements for elliptic interface problems, Appl. Numer. Math., 60 (2010), pp. 19-37. 
[42] I. PERUGIA AND D. SCHOTZAU, On the coupling of local discontinuous Galerkin and conforming finite 

element methods, J. Sci. Comput., 16 (2001), pp. 411^133 (2002). 
[43] M. PLUM AND C. WIENERS, Optimal a priori estimates for interface problems, Numer. Math., 95 (2003), 

pp. 735-759. 

[44] M. PROSI, P. ZUNINO, K. PERKTOLD, AND A. QUARTERONI, Mathematical and numerical models for 
transfer of low-density lipoproteins through the arterial walls: a new methodology for the model set up 
with applications to the study of disturbed lumenalflow, J Biomech, 38 (2005), pp. 903-917. 

[45] A. QUARTERONI, A. VENEZIANI, AND P. ZUNINO, Mathematical and numerical modeling of solute dynam- 
ics in blood flow and arterial walls, SIAM J. Numer. Anal., 39 (2001/02), pp. 1488-151 1 (electronic). 

[46] B. RIVIERE, Discontinuous Galerkin methods for solving elliptic and parabolic equations, vol. 35 of Frontiers 
in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 
2008. Theory and implementation. 

[47] W. RUDIN, Functional analysis, International Series in Pure and Applied Mathematics, McGraw-Hill Inc., 
New York, second ed., 1991. 

[48] C. SCHWAB, p- and hp-finite element methods, Numerical Mathematics and Scientific Computation, The 



DG METHODS FOR MASS TRANSFER THROUGH MEMBRANES 



25 



Clarendon Press Oxford University Press, New York, 1998. Theory and applications in solid and fluid 
mechanics. 

[49] R. STENBERG, Mortaring by a method of J. A. Nitsche, in Computational mechanics (Buenos Aires, 1998), 
Centra Internac. Metodos Numer. Ing., Barcelona, 1998, pp. CD-ROM file. 

[50] T. TEORELL, Transport processes and electrical phenomena in ionic membranes. Prog Biophys Biophys 
Chem, 3 (1953), pp. 305-369. 

[51] E. ZEIDLER, Nonlinear functional analysis and its applications II, Springer, New York, 1990. 

[52] P. ZUNINO, Mathematical and numerical modeling of mass transfer in the vascular system, PhD thesis, Ecole 
Polytechnique Federale de Lausanne, 2002. 



